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Abstract 

The assumption of real-number arithmetic, which is at the basis of conventional geometric algorithms, has been 
seriously challenged in recent years, since digital computers do not exhibit such capability. A geometric predicate 
usually consists of evaluating the sign of some algebraic expression. In most cases, rounded computations yield 
a reliable result, but sometimes rounded arithmetic introduces errors which may invalidate the algorithms. The 
rounded arithmetic may produce an incorrect result only if the exact absolute value of the algebraic expression is 
smaller than some (small) e, which represents the largest error that may arise in the evaluation of the expression. 
The threshold e depends on the structure of the expression and on the adopted computer arithmetic, assuming 
that the input operands are error-free. A pair (arithmetic engine,threshold) is an arithmetic filter. In this paper 
we develop a general technique for assessing the efficacy of an arithmetic filter. The analysis consists of evaluating 
both the threshold and the probability of failure of the filter. To exemplify the approach, under the assumption 
that the input points be chosen randomly in a unit ball or unit cube with uniform density, we analyze the two 
important predicates "which-side" and 'insphere". We show that the probability that the absolute values of the 
corresponding determinants be no larger than some positive value V, with emphasis on small V, is O(V^) for the 
which-side predicate, while for the insphere predicate it is Q{V3) in dimension 1, O^V^) in dimension 2, and 
0(V2 In i) in higher dimensions. Constants are small, and are given in the paper. 

1 Introduction 

The original model of Computational Geometry rests on real-number arithmetic, and under this assumption the 
issue of precision is irrelevant. However, the reality that computer calculations have finite precision has raised an 
increasing awareness of its effect on the quality and even the validity of geometric algorithms conceived within the 
original model, in the sense that algorithm correctness does not automatically translate into program correctness. In 



recent years this issue has been amply debated in the literature (see, e.g., [BKM^95, FV93, Yap97|). In particular, 
it has been observed that while some degree of approximation may be tolerated in geometric constructions, the 
evaluation of predicates ( the "tests " carried out in the execution of programs, - suc h as wh i ch-side , incircle , 



insphere -) must be exact to ensure the structural (topological) correctness of the results [ BMS94 , yap97 , LPT96 1 . 

In principle, error-free predicate evaluation is achievable for error-free input operands, if the latter are treated 
as integers and the arithmetic is carried out with whatever operand length is required to express the intermediate 
results. Such safe approach, however, if adopted in its crudest form, would involve enormous overheads and would be 
nearly impracticable, since the execution time of some operation (such as multiplication) may increase quadratically 
with the length of the representation. 
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As a time-saving alternative to exact arithmetic, it has been customary to resort to rounded (approximate) 
arithmetic (for example, floating-point arithmetic). Such practice can be modeled as follows. Evaluation of a 
predicate P typically involves computing the value ^ of some expression, built using rational operations. The value 
of can be mapped to one of three values: positive, zero, negative ( referred to here as the "sign" of /i), which 
defines the predicate P. Let £ denote an evaluator for P, and let n{£) denote the numerical value computed by £. 
In general, £ is an approximate evaluator of ^, so that its use involves the adoption of a device, called a certifier, 
intended to validate t he corr e ctness of the evaluation. The pair (evaluator, certifier) is what has been refered to as a 



filter in the literature FV9J , MN94]. Typically, the certifier for £ compares |/i(f )| with a fixed threshold e{£) > 0. 
If 1/^(5)! > £{£), then the sign of jJ.{£) is reliable. Otherwise, the certifier is unable to validate the result and we 
have a failure of the filter. In such event recourse to a more powerful filter is in order. This suggests the need to 
develop a family of filters of increasing precision (and complexity), to be used in sequence until failure no longer 
occurs. The last item of this sequence is the exact evaluator, for which the certifier is vacuous (i.e., e{£) — 0). Such 
approach, with an obvious trade-off between efficacy and efficiency, embodies the notion of adaptive precision. 

From a practical standpoint it is therefore very important to gauge the efficacy of very simple filters, that is, 
their probability of success. If it turns out that if a filter has a high probability of success, then recourse to a 
more time-consuming filter (or exact computations) will be a rather rare event. Of course, any such estimate of 
efficacy rests on some arbitrary hypotheses on the a priori probability of problem instances. This is an important 
caveat; however, under reasonable hypotheses (uniform distributions), we submit that the obtained estimates will 
be a significant contribution to the assessment of the validity of such approaches. 

In this paper we analyze two specific predicates, "which-side" and "insphere". Both predicates consist of 
computing the signs of appropriate 5 x S determinants, whose entries are specified with a fixed number of bits. 
Depending upon the adopted evaluation scheme (the choice of an equivalent expression for a given function) and 
upon the precision of the input operands (for example, only a fixed-length prefix of their representation may be 
used in the evaluation), only a prefix of the computed value is reliable. This means that if the absolute value of 
the determinant is above a known threshold e, then its sign is also reliable. 

Our objective is therefore two-fold: 

1. To compute the value of the threshold e for a given determinantal evaluation technique; 

2. To compute the probability that the absolute value of the result of the evaluation does not exceed e, i.e., the 
probability of filter failure. 

We recall that when evaluating the determinant of a (5 x 5 matrix A (a 5-determinant, for short) we are computing 
the signed measure of a hyperparallelepiped defined by the S vectors corresponding to the rows of A. Since each 
of the components of these vectors is an integer in the range [— 2*""^, 2''"^ — 1], a generic vector is (applied to the 
origin and) defined by its free terminus at grid points in a 5-dimensional cube Cs of sidelength 2'' centered at the 
origin. 

Our probability model assumes that all grid points within Cs have identical probability. Our analysis aims at 
estimating (a majorization of) the distribution of the volume of the hyperparallelogram described above. To obtain 
the desired result, we introduce some simplifications consistent with the objective to majorize the probability. 
Specifically, while our assumption is a discrete uniform distribution within Cs, we begin by considering a continuous 
uniform distribution within the ball Bs, the 5-dimensional ball of radius 1. The obtained results are then used 
to bound from above the distribution of the volume for uniform density within the cube Cs (the 5-fold cartesian 
product of interval [—1, 1]), and are finally extended to the target case of uniform discrete distribution. We shall 
recognize that the initial simplification (uniform density in Cs) closely approximates the more realistic situation. 

The paper is organized as follows. We begin with the "which-side" predicate (referred to as "determinant"), 
i.e., in Section ^we carry out the probabilistic analysis in Bs, for S — 1, 2, 3 and arbitrary S (detailed considerations 
of the low-dimensional cases has an obvious pedagogical motivation). In Section |^, we extend the results from the 
continuous ball to the discrete cube. In Section ^ we carry out an analogous analysis for the "insphere" predicate, 
which illustrates the adverse effect of dependencies among the determinant entries. Finally, in Section ^ we evaluate 
the precision of determinant evaluation by recursive expansion, and illustrate the efficacy /efficiency tradeoff. 
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2 Probabilistic analysis within unit ball 



Throughout this section we adopt the following notation: We let xi, X2, ■ . ■ , xg be the coordinates of (5-space,and 
pi,P2, ■ ■ ■ ps be S points in the unit ball of dimension S. We also denote by \pi,p2 ■ . - Psl the absolute value of the 
determinant defined by points pi,p2 ■ ■ - Ps- This quantity, which is the volume of the hyperparallelogram defined 
by the origin and by points pi,p2 ■ ■ .ps, will also be denoted as- 

We begin by examining in some detail the cases of low-dimension determinants. Since the analysis is done in a 
visualizable geometric setting (5 < 3), it is preparatory to the more abstract higher-dimensional cases. 



2.1 1-and 2-determinant 

Obviously, if pi is uniformly distributed between -1 and 1, then 



Prob{\pi\ < R) = R 

Less trivial is the analysis of the two-dimensional case. We will study the probability for |pi,p2| to be smaller 
than a constant A when pi and P2 are distributed uniformly in the unit disk. 

Once pi is chosen (due to the circular symmetry, pi is represented by a single parameter ai , its distance from 
the origin), p2 will yield an area between 02 and 02 -I- da2 if it belongs to one of the two strips of width depicted 
in Figure hi. 

We then have 



Prob{\pi,p2\ < A) = / Prob{\pi,p2\ < A \ ai)p{ai)dai. 
Jo 

Since the density function of ai is p(ai) — i2Uix — 2a\, and the density of 02 conditional on ai is p{a2 \ ai)da 
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i4y 1 - (^) ^ (d) is the density in the unit disk, and 4^1- (^) is the total length of the two strips) we 
have 



Prob{\pi,p2\ < A) 



1 pmin{A,ai ) 



2p{a2\ai)aidaida2 



ai—O J a2 — 
A i p\ , 



• 1a\da\ da2 



-da2 



\/ ai^ — a-i^ — a2 arccos 



02 arccos 02 



/ —da2 f \/l — a: 
8 3 /- y ] 

— A\/ 1 — A^ + — arcsin A — —A^ arccos A 



1 . 1 2 

— arcsm 02 02 arccos 02 

4 2 



One can easily verify that the value of the last expression is 1 for A = 1, which is the maximum attainable value 
for the area of the parallelogram. 
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2.2 3-determinant 

Again we assume that points are uniformly distributed in the unit ball, and compute the probability that the 
volume of the parallelepiped defined by 0,pi,p2, and ps is smaller than a constant V > 0. 

To compute this volume, we begin by considering the parallelogram defined by 0,pi, and p2, evaluate its area, 
and then consider the distance of ps from the plane containing the parallelogram. 

Distance Opi is between ai and ai + dai if pi belongs to a spherical crown of thickness dai and area ATia\. 
Therefore the probability density of ai is p(ai) = -^A-ko^ — 3a? (note that is the density in the unit sphere). 

Once pi has been chosen, the area of the parallelogram defined by 0,pi and p2 is between 02 and 02 +da2 if P2 
belongs to the crown of thickness of a cylinder of radius ^ whose axis contains Opi (see Figure ^) . Therefore 
the distribution of 02 conditional on ai is given by 



p{a2\ai)da2 — 



47r ai V \aiJ ai af y \aiJ 



Finally, once pi and P2 have been chosen, the volume of the parallelepiped is between 03 and 03 + da-^ if p3 



belongs to one of the two spherical slices of width parallel to the plane containing 0,pi, and p2, and at distance 
^ from it (see Figure ^). Therefore the distribution of 03 conditional on a2 is given by 

*-i->--i^--'('-(5)>^^i('-(5)>" 

On the basis of this analysis, we can say 

PTob{\pi,p2,P3\ < V) 

"1 /* Q. 1 /• mi 11 ( V, a 2 ) 

/ / p(ai)p{a2 I ai)p(a3 | a2)daida2das 

I / I / P(«3 I a2)v{a2 I ax)p{ax)da2 ) dai ) da^ 

\^ \J 
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Figure 2: For the analysis of the 3-determinant 
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The latter integral is not elementarily computable. Since the integrand is always positive, so is the integral. To 
neglect it corresponds to majorizing the probability, which is conservative for our analysis. Therefore we write 

27iT 9ir , 3ir , 81 2 / T 27 , 27 

Prob(\pi,p2,P3\<V) < V V'lnV V \/ 1 - arccos V V arcsin V 



If we set V = 1, then the neglected term can be evaluated exactly to ^ — 1 by exchanging the order of 
integration, and we correctly obtain the value 1 for the probability. 

Furthermore, using the inequalities 03 < ai and arcsin ^ < ^, the neglected term can be bounded from above 

¥ lo' la o.S'^d.aidas < ^^V^, so we guarantee the tightness of the approximation of the probability by ^j^V 
when V is small. 

The preceding analysis, in its simplicity, reveals the essential items for the evaluation of the relevant conditional 
probability densities. Specifically, referring concretely to the case 5 = 3, due to the assumption of uniform distribu- 
tion of the points in the unit sphere, the conditional probability density p{ai\ai-i)dai of at given ai-i, i — 2,3, is 
proportional (through the value of the density in the unit sphere) to the volume of some three-dimensional domain. 
The latter is a thin crown (of thickness -^^) of a three-dimensional surface 53,i which is the locus of the points at 

a distance between ^ and ^ -I- from the flat J-i-i spanned by variables xi, . . . ,Xi^i. Surface Ss^i has a 
very simple structure. Let (u, v) be a pair of points realizing the distance °' , with it G J-i-i and v £ S3 i. Point 
V belongs to the boundary of a (4 — i)-dimensional ball of radius ^ of which point u is the center: therefore 
this entire boundary belongs to S^^i ( in the discussion above, this boundary consists of a circle for i = 2 and of 
two points for i — 3 ). Moreover, since u belongs to a flat Ti-\, any translate of v within the unit sphere in a 
flat parallel to Ti-\ also belongs to Sz^i. These translates form the intersection of the unit sphere with a flat at 

. 2 

distance °' from the center of the sphere, and therefore are an (i — l)-dimensional ball of radius 



(in the discussion above this ball consists of a segment for i = 2 and of a disk for i — 3). We conclude that 53,i 
is the cartesian product of the boundary of a (4 — i)-dimensional ball of radius ^ ( the "boundary" term) and 

of an (i — l)-dimensional ball of radius ~' (^T^) i^^^ "domain" term). The expression for the conditional 

probability density consists of four factors: the density within the unit sphere, the measure of the boundary term, 
the measure of the domain term, and the thickness of the crown. These four items (in the order given) are evidenced 
in (|) and (@). 

2.3 Higher-dimensional determinant 

We now extend the preceding analysis to arbitrary dimension 5. If we assume for a one-dimensional volume (a 
distance) the conventional degree of 1, then the volume and the surface of a j-dimensional domain have respective 
degrees j and j — 1. L et Vj (r) and Sj (r) respectively denote volume and surface of a j-dimensional ball of radius r. 



We recaU that[Ber87, 9.12.4.6] 



r ^ ^ f • ( ^ 2V^(^)! . . . , , 
Vi(r) — —j—r tor i even Vi(r) = r tor i odd. 

V- 

The probability density p{ai)dai — Prob{ai < \Opi\ < ai + dai) is obviously given by ^j^dai. Referring 
next to the observations at the end of the preceding subsection, the conditional probability density p{ai \ ai-\)dai 
of ai after pi,P2, ■ • ■ have been chosen (conventionally in the flat described by coordinates X\,X2, ■ ■ ■ , Xi-i) to 

realize the value a^-i, has the following expression: 

p{ai\ai-i)dai = — \-- ■ ss-i+i ( ) ■ Vi^i 
vs(i-) V tSi-l 




Therefore, since Vi{r) = fi(l) • r* and Si{r) = i ■ Vi{l) ■ r' ^ , we have 
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p{ai\ai-i)da. 



-1 



s 



dai 



In the rightmost term above the product of the powers of the a^'s simphfies to , so that 



i-1 



1 _ I I I da. 



The last expression contains a constant depending only on 5, which will be denoted 

Now we can write the probability for the absolute value of the determinant to be no larger than V 

/I i-ai /•as-2 /■min(V,aa_i) * / / / \ \ 

/ ... / / TT Jl- — dai (c) 

/ / I \ \ / I N 

1 /"(il r^d — 2 / /'niin( V,a5 _ ^ ) 11 / \ \ \ I j / \ ^ 



We now observe that the integral 



is trivially bounded by V. Therefore we obtain the following upper bound: 

i-l 

2 ' 

ai 



Prob{\pi,p2...ps\<V)<ksV J ... j Yl^^l- i^-J-j j da. 

Since the latter integral does not depends on V , its value is another constant, which we denote Is and which 
depends only on 5. However, when V — 1, Equation yields Prob{\pi,p2 . . .ps\ < 1) = ksTs+i and, since 
Prob{\pi,p2 . . .psl <!) = !, we get Is = 1777- 
Finally, we conclude 

Prob{\pi,p2...ps\<V) < ksIsV 

< S'-^V^asV 

vs(iy ^ 
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For small S, the bounds are given below. For S = 1, 2, 3 the bounds coincide with the values previously found, that 

is. 



0-1 = 1 
2 



3 



fT2 = — « 2.5 

TT 

= ^«5.3 

213 

-4 = 

224 

3 From continuous ball to discrete cube 
3.1 From continuous ball to continuous cube 

The above calculations have been carried out for points uniformly distributed inside the J-dimensional ball of 
radius 1, referred to as Bs- This assumption may not seem to model the real situation for two reasons: (i) points 

manipulated by computers have discrete rather than continuous coordinates, and (ii) points are more reasonably 
assumed to be uniformly distributed in a cube than in a ball (as in the case when each coordinate is independently 
and uniformly selected). 

We will first show that the previous result relative to the ball Bs induces a similar result for uniform density in 
the unit cube Cs — [—1, 1]*. 

Note that Cs is contained within a (5-dimensional ball of radius \/S, denoted \/5Bs- First, we consider points of 
Cs as points of VSBg and apply a homothety with a factor thereby obtaining 

Prob{\puP2 ■ ■ -Psl <V\pie VSBs) = Prob{\pi,p2 . . .pel < v^"V | pi € Bg) < -^V 

Next, we wish to restrict the points to belong to Cs, i.e., we consider the event pi € VSBs, i ~ 1, . . . ,S, as the 
union of the event pi € Cs,i = 1, ... ,5, and its negation. The probability of this event is clearly ( — ^-—i=j] , so 
that 

Pro6(|pi,p2 . . ■ps\ <V\pi^ VSBs) > i ^— ^ ) Prob{\pup2 . ..psl <V\pieCs) 

\vs{l)Vd'J 

which gives us an upper bound to the probablity in question. Specifically 

vsWVS \ as 5vsil)vs-iilfVS a , 
2^ 1 ^ " 

Practically, for small values of 5 we get: 

V-i = 1 

1p2 = 



8 



128 

V4 = ^71-' « 380 

/ 9765625 12 _ „„„„„ 

= 402653184 " "^3000 

19683 15 ^., 
= 125000" 

The previous computation can probably be generalized to other kinds of domains provided that the ratio between 
the volumes of the inscribed ball and of the circumscribing ball is bounded. 



3.2 From continuous cube to discrete cube 

In this section we shall discuss why the obtained results for continuous density of points are still useful for discrete 
probability, i.e., when points belong to a regular grid of ^ points inside Cs- 

Notice that the preceding results are clearly incorrect for discroto probabilities. In fact, they prescribe: 
Prob{\pi,p2 ■ ■ - psl = 0) = 0, which is false for discrete probabilities. For example, in two dimensions, when pi 
is chosen, p2 coincides with pi or with the origin with probability 2if, and yet the determinant value is 0. For 
a less trivial case (when pi,P2, and the origin are distinct), the event pi = {xi,X2) and p2 — (2a::i, 22:2), with 
xi < and yi < ^, has probability > while the determinant is still 0. However we can still prove that 
Prob{\pi . . - ps] <V)= 0{V) when 77 is smaller than V. 

If pi .. .pa is a set of points in Cs whose determinant is not larger than V, we will map it to a the nearest set of 
grid points p'l . . .pg whose determinant is not too large. More precisely, ii V = (pi, . . . ,ps), V = (p'l, . . . ,Pa) and 
V = 'P' -V = {di,...,di),vfe have 

\p'\ = \P + V\ = \P\+ ^ \VVi\ 

where / is a nonempty subset of {1, 2 ... ,5} and \'P'Di\ is the determinant obtained by replacing, for each i £ I, pi 
with di in \V\. The above result follows from the multilinearity of the determinant. If the cardinality of I is j, we 
can bound from above the absolute value of \VDi\ by the product of the norms of its vector. Since \\pi\\ < \f6 and 

11*11 < we have \VVi\ <V~& ^.By grouping the (^) terms with identical value of j we get: 



\V'\-\V\ 



<Vd' 



Referring now to the ^^-dimensional space whose points are the sets pi . . .ps, the determinant \V'\ is no larger 

than V if all the points V in the hypervoxel (in 5^ dimensions) have determinant IV] no larger than V + SVS ^. 
Since, clearly, a random point V for the continuous distribution can be in any voxel with the same probability, we 
conclude that 



Probl |pi...p5| < V 



discrete 
distribution in 



Cs < Proh [ \pi...p5 



continuous 
distribution in 



Cs 



< V5.(V + ,5V^'|) 
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4 The insphere test 



In the preceding analysis of the "which-side" predicate, the points defining the hyperparallelogram were assumed 
to be independent and equally distributed. In this section we consider a case for which there exist dependencies 
among the coordinates of the points: the "insphere" predicate. This predicate, referred to as (5-insphere for short, 
tests whether in S dimensions the origin lies inside the hypersphere defined by (S + 1) other arbitrary points 
Pi = {xi,i,Xi^2, . • . , Xi^s), i = 1,2, . . . ,5 + 1. 

It is well known that the 5-dimensional insphere test is embodied in the sign of the determinant 



Xll 
X21 



Xl2 
X22 



Xs+2,1 Xs+2,2 



xll + T?2 + . . . + x\s 
xll + XI2 + ■■■+ xlg 



Xs+2.1 + Xs+2,2 + . ■ . + Xs^2,S 



Without loss of generality, one of these points {ps+2) can be chosen as the origin O, so that the above determinant 
becomes 



Xll 

X21 



X12 
X22 



2:5+1,1 xs+1^2 



■^11 
xll 



2 I 2 

Xs+1.1 + 3:^+1,2 



" xis 

' X2S 



+ Xs+l,S 



These 5 + 1 points are assumed to be evenly distributed in the unit cube Cs 



4.1 1-insphere test 

We can model the problem as follows. The origin O and u define the 1-dimensional sphere; v is the query point. 
Parameters u and v are independent and uniformly distributed in [—1, 1], and we wish to evaluate the probability 
of the following event: 



lAil < A 



with Ai = 



We can view the above determinant as defined by two points {u,v) and in the u,v plane. Clearly 

the choice of point (u, v) completely determines the determinant value. Point {u, v) is uniformly distributed in 
the square [—1, 1] x [—1, 1]. Since uv{v — u) is the determinant value, the determinant is null on the three lines 
u = 0, V — and u — v; its values are symmetric with respect to the line it — —v and antisymmetric with 
respect to u = v. Therefore for any value of < yl < 2 it is sufficient to evaluate the probability in the quadrant 
—u < V < u,0 < u < 1 (fully shaded quadrant in Figure ^) and multiply it by 4. In the upper semiquadrant (where 
the determinant is negative) the contour lines have equations: 



.^1 
2 2 



.A 



The two curves for fixed A join with a common vertical tangent at the point [n, ^) where /i = (4 A) 3 (notice that 
such curves exists only for A < i). In the lower semiquadrant (where the determinant is positive) contour lines 
have equations: 

u 1 n, 7a 

V — 1 

2 2 



+4- 



and 



and intersect the line v = —u at point (u,—!/), where u = i(4y4)3 



The probability of the event described 



above is given by the heavily shaded area in the Figure ^ (in fact this area should be multiplied by 4 since there 
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are four quadrants, and normalized, dividing by 4, since 4 is the area of the square). The area is also 1 minus the 
area of the lightly shaded region. The latter area, which we want to bound from below is given by: 



1 - 4:4rdu 



We now observe that 



and (for A < ^] 



so that 



I > 



l-44>l-44 



udu — A 
1 



udu — 4A 



-du 



2-T + ^-- 



1 M . . 4^1 
— + 44 

2 2 



If we now use the fact the fact that 2u — fj, = (4^4) a , then we obtain: 



Prob(|Ai| < A)< 



A3 - < 5.355A3 



A direct numerical calculation gives Pr-o6(|Ai| < \) = 0.7, while the value of the above bounding expression 
for the same value of A is 0.85 ( an excellent agreement considering the rather high value of A). For A > i, the 
lightly-shaded region in the upper semi-quadrant disappears yielding 

Prob(|Ai| <A)< 3v^A* - 44 < 3.78A3 

but, considering the high range of A, this expression has little practical interest. 

Alternatively, we may consider the following formulation. Rather than lifting point v to the parabola y = 
v^u — u^v (for fixed positive it), we lift it to the parabola y = — uv = F{v) (which intersects the i7-axis at and 
u) so that 

Pro6(jAj <A) = Prob{\F(v)\ < -) 



This formulation will be useful when considering the multidimensional case in Section 4.2.2 



4.2 Higher dimensional insphere test 

By elementary column operations the determinant A^ can be transformed into one where the last column has zero 
entries except in the last row. Specifically, A^ = defines the circumsphere 5* of the points 0,pi, . . . ,ps, whose 
center is the point (^, . . . , ^). Subtracting column i times Ci {i = 1, . . . , 5) from the last column, we obtain: 



As = 



111 

X21 



X12 
X22 



XU 
X2S 



Xsi XS2 

a;5+i,i xs+1,2 



xss 

Xs+1,5 W 
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Figure 3: The (w, u)-region defining (|A| < A) in the 1-insphere test 



The determinant of the intersection of the first S rows and columns of the above matrix gives the signed volume v 
of the hyperparallelogram defined by the first 5 points and W has the dimension of the square of a length and is 
the value of the function 

^(2:1, . . . , xs) = xi + xl + . . . + xj — cixi — . . . cgxs (d) 



evaluated at point ps+i- We let = i(cf + ... + c|) denote the square of the radius of sphere S. Notice 
also that xs+i — F{xi, . . . , xs) is the equation of a hyperparaboloid 7i in (5 + l)-dimensional space, so that 
F{xs+i,i, . ■ ■ , xs+i,s) = F{ps+i) is the signed height of the point obtained by lifting ps+i from the hyperplane 
2:^+1 — (to which it belongs) to Ti. Notice that hyperplane xs+i — contains the hypersphere S. 

We now wish to bound from above the probability of the event \F{ps+i)\ < V , for some constant V . 

Assuming as usual constant density, this probability is the volume of a sphere S" of radius u" = \/V + v? minus 
the volume of a concentric sphere 5*' of radius u = \/ —V + u^, whenever the latter is defined. (These spheres are 
intersected with Cs and normalized by its volume.) In general dimension, since the above radii depend upon the 
points pi,p2, ■ ■ ■ ,ps, the evaluations of the volumes of the above spheres is problematic. However, we shall show 
that an interesting simplification occurs for 5 = 2. 

Below we shall use the following bounding technique. Let As be expressed as the product of two continuous 
random variables 

= ah 

Given a constant a we have: 

V 

Prob{ab<V) < Prob{a < a or b < —) 



that is: 



V V 
< Prob{a <a) + Prob{b < — ) - Prob{a < « and 6 < — ) 

Q a 



Prob{ab < V) < Prob{a < a) + Prob{b < — ) 



(e) 
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4.2.1 2-insphere test 

In this case the volumes of S" and S' are respectively Tr{V + u^) and ■k(—V + u^), so that, if S' is defined, then 
their difference becomes 7r2V, otherwise (i.e., when < V) tt{V + u^) < tt2V. After normalization (since 4 is the 
measure of the unit square) we have 

Prob{\Fip3)\ <V)<^. 
Given that IA2I — \piP2\F{p3) inequality (^ becomes 

Prob(\A2\<V) < Prob(\pip2\ < a) + Prob(\F(p3)\ <—) 

a 

^ , ttV 
< ip2a+——. 

2a 

From Section 3 we know that -02 = ti". Selecting for a the critical value a — we have 

Prob(| A2I <V) = nV2V ^ 4.44W 

4.2.2 5-insphere test (5 > 2) 

As for the two-dimensional case, we shall prove that, in general, the probability that \F{ps+i)\ is no larger than W 
is sufficiently small. 

Again, inequality (^) becomes: 

W 

Prob(\As\ <W)< Prob{\pi...ps\ < —)+ Prob{\F{ps+i)\ < a) (f) 

Q 

We know, from the results of Section ^, that 

W W 

Prob{\pi . . .ps] < — ) < ips — 
a a 

Therefore, there remains to bound from above Prob{\F{ps+i) \ < a) . Recall that is the distance from 

the plane xs+i = of the point ps+i lifted to the hyperparaboloid Ti, of equation F{p) — 0, which intersects 
xs^i = in a 5-dimensional sphere passing by the origin with with center q = (^, . . . , ^). We define point 
such that it lies on a line I passing by O and q and such that length(p|^j^, q) — length(p^+i, q) and among the two 
possible choices, we select the one closest to O (see Figure ^). It is immediate that |F(pi+i)| = 1-^(^1+1)1 and that 
length(Op|+i) < Vs. 



Thus, our problem is reduced to a one-dimensional instance, closely related to the one we studied in Section 4.1 
(here u — 21ength(Og) and v is the signed value of length(Op|_|_i)). There are however, some significant differences. 
There, we were evaluating the probability of the event |'!i-F(w)| < V , and variables u and v were uniformly distributed 
in [—1, 1]. Here, on the other hand, we wish to evaluate the probability of the event |-F(«)| < a, "radius" u varies 
between and 00, and "distance" v varies between and Moreover, the densities pi{u) and P2(v) of u and 
V respectively are not constant; however, as we shall see, they are bounded by constants qi and 52. 

Next, we observe that 

Prob{\F{ps+i)\ < a) ^ pi{u)p2{v)Prob{\F{ps+i)\ < a \ u, v)dvdu 

■J It— -J v=—\/5 

Pi{u)p2{v)Prob(\F{ps+i)\ < a \ u,v)dvdu 



Jv=0 

2 



Since the function F{ps+i) has the expression v — uv, the right-hand-side of the above equation is just the integral 
of pi{u)p2{v) on the domain( shown in Figure ^ bounded by the curves v'^ — uv = ±a. This integration domain 
can be split into four subdomains A,B,C,D, as shown in Figure ^. This split depends on a parameter f3 < -Ja which 
will be chosen later. In detail we have: 



13 



Figure 4: Definition of p|_|_x, u and v, using a 2-dimensional instantiation 



A 



Pi{u)p2{v)dvdu < qiq2 I — du < (jri(j'2 ( In 5 + 2 In — 1 a 



Domain A is defined by /3 < u < yj\5) and v — ^ <w<i;+^. pi{u)p2{v) < qiq^ and thus 

— fJ.Tj < «i (7o I In rl 4- 9 In 

P. 

Domain B is a rectangle defined by < w < /3 and (3— ^<u<l3+^. pi{u)p2{v) < qiq2 and thus 

J J pi{u)p2{v)dvdu < <ji(72/3^ < 2qiq2a 

Domain C is defined hy < v < P and /3 + ^ ^ u < v + ^. P2{v) < q2 and since on any 

domain, we get 



Pi{u)p2{v)dvdu < q2l3 



C 
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Figure 5: Integration domain 
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Figure 6: For the analysis of pi{u) 



• Domain D is analogous to domain C, 




Pi{u)p2{v)dvdu < q2l3 



These results yield: 

Prob{\F{ps+i)\ < a) = 2giQ2a ^lnj + 21ni +2^ +4q2(3. (g) 
To obtain values for gi and g2, we observe: 

1. With reference to pi{u), consider the (5— l)-dimensional sphere passing by 0,p\. . .ps-i in the flat fl defined by 
these points (see Figure ^ for an illustration) , and let t be its center. Consider the family of (5-dimensional 
spheres passing through the origin and whose centers project to t in fl. Point ps will determine a (5-dimensional 
sphere of radius between u and m + dit if it lies between a pair of spheres of (or its symmetric pair with 
respect to 11) of radii u and it + du. The volume p\{u)du of this region (normalized by the volume 2^ of Cs) 
is < -^Asdu, where As is the maximum of the measure of the surface of a 5-dimensional sphere passing by 
the origin and intersected with Cs. A trivial upper bound on As is the measure 252^'^ of the surface of Cs- 
Therefore 

p,[u) < ^25^51 

2. With reference to P2{v), we observe that once ps is fixed, the choice of ps+i will produce a value between v 
and V + dv if ps+i belongs to the shaded region in Figure vA The volume of this region is dv times the surface 
of the sphere of radius v inside Cs ■ This surface is clearly less than the surface of Cs , so that we get 

P2{v) < ^.252^-1 =5^q2 
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Figure 7: For the analysis oi p2{v) 



Substituting gi = 25 and q2 ~ 5 into (jg) and choosing /3 = ^ < ^/a we obtain: 

Prob(\F(ps+i)\ <a)< 85^aln-+ (4(5Mn 5 + 85^ In 2 + 85^ + 5) a = r^alni + Osa 

a ^ 'a 

(Small values of ts and are rs = 72, 6*3 = 164.5, T4 = 128, 6*4 = 309.4, rg = 200, 65 = 504.6, re = 288 and 
06 = 751.6.) Inequality (|) becomes: 

W 1 

Prob(\As\ < W) < ips — +riQln- + Osa 
a a 

so that, setting a = \J^^^^ 

Proh{\As\ <W)< y'Mn + es)^ (In 7^ - In 1±±1± + 1 ) 

\ W -ips J 



<t>3 ~ 70 


X3 


= -100 


04 ~ 408 


X4 


= 350 


05 « 3970 


X5 


= 18000 


6 ^ 68500 


X6 


= 640000 
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4.3 Discrete distribution 



As for the case of the which-side test, also for the insphere test we can map each point to the nearest gri d po int, 
and evaluate the ensueing effect on the determinant to be computed. Generalizing the notation of Section 3.2, we 

have that the norm j|p]"|| of the lifted point is bounded by \^T+W and that < V5^2^ + 2| < ^v^r; + rj'^. 
Therefore, grouping the errors term, we get 

— -. .s^i 

-s+i 



|7?'t| _ 


-\v'\ 


< 









v2 



■(5 + 1)^5 + 52 



■5+1 



y/S + 5- 



and considering rj small 



5 The efficacy of arithmetic filters 

To complete the analysis, in this section we wish to assess, under the given statistical assumptions, the probability 
of failure of a given arithmetic filter for determinant sign evaluation (i.e., the probability that the filter is unable 
to certify the correctness of the computed sign). This probability is a quantitative measure of the efficacy of the 
filter. If |pi . . > es, then the result of the computation is reliable, and so is its sign. Plugging eg in place of V 
in the expression above, we obtain the condition: 



Prob |pi . . .p^l < £s 



discrete 
distribution 



where we have used the result (Section 3) that Pro6(|pi . . .ps\ < V) < ipsV. The parameter ps introduced here is 
therefore the sought measure of filter efficacy. 

To exemplify this approach, we shall compute ps for the evaluation of determinants, for the case where the 
coordinates of the points are floating-point numbers in the interval [—1,1], the computations are carried out using 
floating-point arithmetic with b bits of mantissa, and the determinant is evaluated by standard expansion with 
respect to one of its columns (recursive evaluation). 

To this end, it is necessary to compute the parameter eg- We introduce the following notation: 

• £ [M, m] denotes the set of numbers whose absolute value is bounded by M and whose error is bounded by 
m. Original entries belong to £[1,0]. 

• M denotes 2^'°'5A/1^ 

With this notation, if xi £ £[Mi,mi] and X2 G £[M2,m2], then we have: 

xi+X2e £[Mi + M2,2-^~^ ■ Ml -I- M2 + mi -f 7712] 

xi • a;2 G £{Mi ■ M2,2''''^ ■ Mi ■ M2 + rm ■ M2 + m2 ■ Mi] 

These rules express the mechanics of 6-bit mantissa normalizing floating-point operations with round-off (round-off 
is done to the nearest, thus the error done is half of the value of the last bit). After transforming variable y to an 
£[M,m] pair, we shall express the above rules as an arithmetics on such pairs, as follows: 

1. £[Mi,mi]+£[M2,m2] = £[Mi + M2,2-''-^ ■ Mi+M2+mi+ ma] 

2. £:[Mi,mi] ■ £:[M2,m2] = £[Mi ■ M2,2-''-^ ■ Mi ■ M2 + rm ■ M2 + ma ■ Mi] 
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The results given below are obtained in an Appendix to this paper by the mechanical application (with a few 
noted exceptions) of the above two rules to the recursive evaluation of a determinant. Using a &-bit mantissa, the 
results are: 



• 


£2 


= 2-2-'' 


• £5 


= 576 • 2-^ 


• £8 = 226624 • 2-'' 


• 


£3 


= 13 • 2-^ 


• £6 


= 3672 • 2-'' 




• 


£4 


= 76 ■ 2-^ 


• £7 


= 27304 • 2-'' 





For example, using the IEEE norm on 6 = 53 bits, we can therefore estimate the corresponding probablility of 
failure. The pertinent values of £5, 5\/6 ^, and ps are displayed below: 



• 


£2 


= 2.2- 


10-16 


• 




= 1.6 • 


10" 


16 


• 


P2 


= 1.2 • 


10" 


15 


• 


£3 


= 1.4. 


10-is 


• 




= 8.7. 


10" 


■16 


• 


P3 


= 4.8. 


10" 


14 


• 


£4 


= 8.4. 


10-is 


• 




= 7.1 • 


10" 


■15 


• 


P4 


= 5.9. 


10" 


12 


• 


£6 


= 5.7. 


10-14 


• 




= 7.8. 


10" 


■14 


• 


P5 


= 3.0. 


10" 


9 


• 


£6 


= 8.3. 


10-" 


• 




= 1.1 • 


10" 


■12 


• 


P6 


= 8.7 • 


10" 
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Finally, it is important to evaluate the efficiency of the examined filter. We observe here that a recursive 
evaluation uses 6\ operations, but the cost can be reduced using dynamic programming. In fact, the recursive 
evaluation involves minors of dimension i. In turn each such minor involves (2i — 1) arithmetic operations {i 
multiplication and i — 1 additions), whose operands are either minors of smaller dimension or original coefficients 
(for i = 2, where the recursion stops). Thus, the total number of operations is: 

= E (0 (^C-^ - 3) - 1) = (5 - 1)(2' - 1) ~ 52'. 
^ ^ 

For 5 < 7 we obtain n = 0, r2 = 3,r3 = 14, n = 45, rs = 124, re = 315,r7 = 762 and rg = 1785. 

It must be pointed out that, for the same function considered above (determinant value), one may choose 
alternative evaluation schemes (corresponding to different expressions of the given function, such as other expansion 
rules, Gaussian elimination, etc.) and/or different arithmetic engines. Each such choice would embody a filter, 
whose efficacy and efficiency can be assessed with the outlined method. 

6 Summary of results and conclusions 

In this paper we have developed a general approach to the assessment of the efficacy of arithmetic filters, under 
some reasonable probability assumptions. As an important example, we have considered the efficacy of filters for the 
evaluation of signs of determinants, both for a case where all entries are independent (the "which-side" predicate) 
and for a case where dependencies exist (the "insphere" predicate). This analysis, in general, consists of two parts. 

The first part aims at computing the threshold for certification by the filter (i.e., the meiximum error which can 
be generated by the evaluation process) , and does not rest on any assumption about the distribution of the data. 
As an example, we have carried out this threshold analysis for the so-called recursive evaluation procedure, which 
computes a determinant by expanding it with respect to one of its columns. 

The second part, which is considerably more subtle, aims at establishing the probability of failure of the filter, 
i.e., the probability that the result of the computation falls below the threshold. This analysis rests on a priori 
assumptions on the distribution of the input data, which we have taken as uniform within their representation range, 
and has been carried out for the two important geometric tests mentioned above. With the notations introduced 
in the preceding sections, the results are summarized below. 
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Prob{\pi,p2 ■ . -Psl < V \ Pi continuous in Bg) < crgV 
Prob{\pi,p2 • • - Pal < V I Pi continuous in Cg) < ^gV 
Prob{\pi,p2 ■ ■ -Psl <V \pi e Jj-gridinQ) < ipg-^V + agT]) 

Prob{\Ai \ <W \ pi continuous in Ci) < 5.36VKt 

Pro6(| A j| < W I Pi continuous in Cj) < ^gVW In— +Xd^^ 

W 

Prob{\Ai\ < W \ Pi e v-grid in Cg) < 4>6 \/W + ft0?ln / +Xi \/W + PsVv 

Below we repeat for synoptic convenience the definitions of the relevant constants and a tabulation of their 
values for small S. 

.vg.^ilf 
erg = 



i>6 



vg{iy-^ 

5vg{l)vg-,{lY^^'-'^ 



(t>g = ^ (252 In <5 + 4J2 In 2 + 8<52 + 

/ 252ln5 + 452ln2 + 852 + | 
X5 = <;/)5 1 - In 



ag = 



sVd' 



■s 



136 

The values for small 6 are 



5 


erg 






XS 


ag 


Pg 


1 


1 


1 






0.5 


2 


2 


2.5 


3.2 




4.4 


2 


18 


3 


5.3 


21 


70 


-100 


7.8 


200 


4 


10 


380 


408 


350 


32 


2800 


5 


19 


23000 


3970 


18000 


140 


47000 


6 


35 


4.5 • 10** 


68500 


640000 


648 


900000 



The above values are tight only for the constant a. Due to data interdependencies, the analysis of the in- 
sphere predicate is considerably more involved than that of the which-side predicate, and the adverse effect of the 
dependencies is manifest in the larger values of the probability of failure. 

Such analysis is particularly valuable for estimating the time required to test the determinant sign. Under the 
given probability assumptions, we may conclude that for small dimension (< 6) straightforward floating-point fllters 
(i.e., floating-point evaluators) axe extraordinarily effective. 
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7 Appendix. Error evaluation 

Di represents a generic determinant of dimension i, and x a generic original coordinate. It must be pointed out 
that the recursive evaluation technique yields an upper bound of 5\ for Ds because dependencies among data are 
not exploited. However,smaller values of such bound are known: Since the determinant value is bounded by the 
product of the norms of its 6 components, y/S is an upper bound on Ds, which is attained when 5 is a power of 
two (Hadamard matrices). For small values of S the following bounds have been obtained ( some by exhaustive 
calculation): Di < 1, D2 < 2, D3 < 4, L»4 < 16, D5 < 48, De < 160, D7 < 576 and Ds < 4096. Whenever 
applicable, we shall use these results below to obtain tighter estimates. These estimates have the form 

Ds € £lGs,es] 

where es is the object of the analysis and Gs, the largest value attainable by Ds, is only needed to carry out the 
analysis. 

1. Di G £[1,0], by definition. 

2. D2 € S[2, 2-''+^]. In fa.ct D2 = X ■ Di + X ■ Di, X & £[1, 0] and Di e £[1, 0]. Therefore 

5(1,0] • 5(1,0] +£[1,0] -ffl.O] = £[1,2-''-'^] +£[1,2-''-'^] =£[2,2-''+'^] 

3. D3 G 5[4, 13 ■ 2^'']. In fact, D:i = {x ■ D2 + x ■ D2) + x ■ D2. Therefore: 

D3 e {£[l,0]-£[2,2-''+'']+£[l,0]-£[2,2-''+''^])+£[l,0\-£[2,2-''+'-] 
= £[4, 2-''-i • 4 + 6 • 2"'] + £[2, 3 • 2"''] 
= £[4, 2-''-^ • 4 -h 11 • 2"''] = £[4, 13 • 2"''] 

where we have used the fact D4 < 4. 

4. L>4 € £[16, 19 • 2-*'+^]. The result is obtained by applying the previous rules and x € £[1,0], D3 € £[4, 13 • 2"''] 
to the following evaluation scheme: 

D4={X-D3+X-D3) + {X-D3+X- D3). 
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5. Ds G f [48, 129-2 "'"^]. The result is obtained by applying the previous rules and x € f[l,0], D4, £ 
f [16, 19 • 2"''"'"^] to the following evaluation scheme: 

= {{x ■ D4, + X ■ Da) + X ■ Da) + (x ■ Da + x ■ Da). 

making also use of the fact Ds < 48. 

6. Dq £ f [160, 459 • 2"''+^]. The result is obtained by applying the previous rules and x € f[l,0], D5 € 
f [48, 125 • 2"''"'"^] to the following evaluation scheme: 

De = {{{x ■D5+x-D5) + {x-D5+x- D5)) + {x ■ D5 + x ■ D5)) 

making also use of the fact De < 160. 

7. Dr e £[576,3413 • 2-''+3]. The result is obtained by applying the previous rules and x € £[1,0], De € 
f [160,459 ■ 2"''"'"^] to the following evaluation scheme: 

D7 = {{{x ■D6 + x-D6) + {x-D6+x- De)) + {{x ■ De + x ■ De) + x ■ De)) 

making also use of the fact Dr < 576. 

8. Da e £[4096,3541-2-''+*']. A gain, the result is obtained by applying the previous rules and x € f[l,0], 
D7 € £[576,3413 • 2"''+^] to the following evaluation scheme: 

Ds = {{{x ■ D7 + X ■ D7) + (x ■ D7 + a; • D7)) + {{x ■ D7 + x ■ D7) + {x ■ D7 + x ■ D7))) 

making also use of the fact Dg < 4096. 

REMARK. If the original values do not belongs to £[1, 0] but to £[1, e] and e is small enough, then the final error on 
Di must be increased by i\c. c small enough means that the error on M docs not affect the value M, in particular 
if all the manipulated sets are of the form £[M,m] with M = M + m, then adopting a first-order approximation is 
legitimate. 
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